Simulation of open quantum systems 
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We present an approach for the semiclassical treatment of open quantum systems. An expansion 
into localized states allows restriction of a simulation to a fraction of the environment that is located 
within a predefined vicinity of the system. Adding and dropping environmental particles during the 
simulation yields an effective reduction of the size of the system that is being treated. 
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Any rigorous treatment of open system dynamics is 
impeded by the de facto impossibility of simulating the 
dynamics of a system comprised of the system of interest 
and the typically macro-, or at least mesoscopic envi- 
ronment. Numerous prior approaches end up with ef- 
fective equations of motion for the system itself, where 
the dynamics of the environment is eliminated, such as 
a Lindblad equation [TJ [5] that describes the influence 
of an environment onto a system. Such an effective the- 
ory, however, is always subject to approximations, such 
as a Markov approximation, and/or taking the system- 
bath interaction, and/or the environment itself to be har- 
monic. 

In this paper we describe an approach to open system 
dynamics based on the propagation of Gaussian states. 
That is, an initial state is expanded in terms of Gaussian 
states, and each of these components is then propagated 
individually. Under the evolution of a general Hamilto- 
nian an initial Gaussian state will not remain Gaussian. 
However, for well localized wave packets the Hamilto- 
nian can be approximated to be locally harmonic, so 
that the Gaussian character is preserved [21 IH IS] ■ Such 
techniques have proven to accurately describe dynami- 
cal properties of various closed systems, thermodynamic 
properties [BJ |7] , and also open quantum systems [5J |3] • 

We make use of the expansion into localized states 
which permits effectively reducing the size of the environ- 
ment to be considered. The underlying idea is to expand 
the state of the entire system (that is system and environ- 
ment) in terms of Gaussian states, such that the actual 
dynamics can be reduced to those environment compo- 
nents that are situated in close proximity to the system 
particle. This allows elimination of all other environmen- 
tal degrees of freedom from the problem. On the other 
hand, there will be environment particles that approach 
and enter this vicinity, and their degrees of freedom will 
need to be incorporated in the dynamics as sketched in 
Fig. [I] 

Any Gaussian state is completely characterized in 
terms of the expectation values of all coordinates and 
momenta x = [r-y , p\ , r 2 , P2 , • ■ ■ , r n , p n ] = Tr xg , and the 
corresponding covariances 




FIG. 1: Schematic representation of the system particle (large 
wave packet) with a Gaussian wave function interacting with 
its environment (small wave packets). The dynamics that 
is relevant for the system can be restricted to the vicinity 
around the particle sketched by the circle. The environment 
particle that leaves the vicinity to the left can be traced over 
since its impact on the dynamics of the system particle is 
negligible. On the other hand, there is another environment 
particle approaching the system particle from the right, and 
needs to be taken into account. 



However, an explicit parametrization of a Gaussian den- 
sity matrix g in terms of those parameters is a rather 
lengthy expression. The Weyl symbol of the density ma- 
trix, i.e. the Wigner function 



W(r,p) = J d n q(f- ±\g\x+ !>eK« (2 ) 

is more convenient to use, since it has a significantly sim- 
pler parametrization: 
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The evolution of a Gaussian quantum state due to a 
quadratic Hamiltonian gives rise to Newton's equations 
of motion x = SWH with the symplectic matrix S, and 
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where TL contains the second derivatives of the Hamil- 
tonian H with respect to the coordinates and momenta, 

Eq. Q describes the unitary evolution of the system 
and its surrounding environment . Any dissipative nature 
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FIG. 2: The Wigner function for the superposition of two 
Gaussian states \a) + \/3) is comprised of four parts. The two 
'diagonal' parts corresponding to |a), and \/3) respectively, 
are proper Wigner functions on their own, and their evolu- 
tion is described by Eq. Q. The Weyl symbol for the two 
'off-diagonal' parts corresponding to \a){(3\, and |/3}(q| are in 
general complex, but they can be propagated similary to a 
regular Wigner function. However, their dynamics is not gov- 
erned by the Hamiltonian around the phase space position, 
at which it is located, but around the positions of both |a), 
and \/3). 

of the the dynamics is due to tracing over environmental 
degrees of freedom. Obviously, the expectation value of 
any observable A on particles that are not being traced 
over are unaffected by the partial tracing 

TvAp = TvA{Tv p p) . (5) 

Therefore, tracing over some degrees of freedom is equiv- 
alent to simply dropping all components of x and S corre- 
sponding to the degrees of freedom being traced over. In 
turn, adding a particle corresponds to extending x and £ 
by the respective quantities, with vanishing correlations 
between the added particle, and the residual system. 

So far, we have been assuming the initial state be a 
spatially narrow Gaussian wave packet. This condition, 
however, is often not satisfied, and for wider wave pack- 
ets the harmonic approximation breaks down. Neverthe- 
less, any quantum state g (pure or mixed) can be de- 
composed into an incoherent mixture of coherent states 
g = J dp{a) P a {g) |a)(ck|, with the P-function P a (g). 
Since, this holds also for pure states, P a cannot always 
be a probability distribution, i.e. it can adopt negative 
values. The central merit of this representation is that 
different initial states \a) can be propagated individu- 
ally, and the overall final state is then the incoherent sum 
over the individually propagated states. In the presently 
discussed framework of open system, this implies that 
partial traces can be taken at any instance. The big 
disadvantage of such a representation is the often wild 
behavior of P a . In particular for non-classical states P a 
is rapidly oscillating, and close to singular, what severely 



limits its usefulness for practical purposes. 

Often it is significantly easier to expand the state \^f) 
into a coherent sum of Gaussian states \^f) ~ J2 a ^a\ a )- 
Due to the linearity of the Schrodinger equation, one can 
of course propagate each initial term \cti) individually, 
and, then reconstruct the final state 

W t |*) = Y,*cUt\a) * 5^*aW t a |a) , (6) 

a a 

where lAf = T exp(—i/h J dtH a ) is the approximate 
propagator, and H a is the locally quadratic approxima- 
tion of H expanded around x Q . 

In order to take partial traces, however, one needs to 
consider the corresponding density matrix: 

g t = Ut\a)({3\U} c^^Mp K\a){f3\uf ■ 

a/3 a/3 

(7) 

Here, U^' 13 is the unitary transformation generated by the 
Hamiltonian H expanded around the classical positions 
of the states |ev) and \0) respectively. Thus, in the follow- 
ing, we will consider the evolution of the individual oper- 
ators p a p resulting from |a)(/3|. Since those operators are 
not necessarily normalized, all expectation values will be 
defined including normalization (A) p — TrAp/Ttp. Do- 
ing so, one obtains the equations of motion 

f)T ih 2? 

— = 2(SH+5 - - -jSHS - j^H-Z , (8) 

with H + — l/2(7i a + Hp), and H- = H a —Hp; and with 
H a , Hp being the second order expansion coefficients of 
H taken along the phase space positions (a|x|a), (/3|x|/3) 
as depicted in Fig. [2j The equations of motion for the 
positions and momenta x are a bit more lengthy, and 
they are most conveniently characterized by the relation 

1 i 
x = -(x Q + x^) - -S5(x Q - ■x.p) , (9) 

and x a and evolving according to Newton's equation 
of motion. 

The Weyl symbol corresponding to those complex 
phase space coordinates and uncertainties is given by 
Eq. ([3| up to the additional factor exp(?7) with 

1 i 
V= ^(x Q -x )5S5(x Q -x^)-— (p Q +p /3 )(f a -r / 3) , 

(10) 

and a phase (p whose evolution is given by (p = C a — Cp — 
iTr£7i_, where £j, (i = 1,2) is the Lagrange function 
with variables ?i , and pi . 

So far we were concerned mainly with the unitary dy- 
namics of the system. However, some care has to be 
taken while tracing over an environmental particle. As 
mentioned before, the dynamical variables x, and S re- 
main unchanged under partial tracing. However, x Q , and 
are propagated rather than x, and for a general state 
any entry of x depends on all the entries of x a , and x^ 
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FIG. 3: The coherence of the superposition of |q) and |/3) gets 
reduced in time due to interactions with an environment gas. 
Here, this coherence is characterized in terms of the Hilbert- 
Schmidt norm of the operator that the initial operator |q)(/3| 
evolves into, and it is plotted as a function of time for 4 dif- 
ferent initial states, where the relative phase space separation 
5x between \a) and \/3) is chosen to be 10, 20, 30, and 40 os- 
cillator length. The coherence decays exponentially, and the 
decays is faster for larger initial phase space separations. 



as shown in Eq. Therefore, one may not simply drop 
the entries of x„, and x^ that belong to a particle that is 
traced over, but one needs to consider their contributions 
to x first. More explicitly, after a partial trace has been 
performed, x is not given by Eq. ^ anymore, but there 
is the additional term i/HY.SP(x a — x^) where P is the 
projector onto the the space spanned by the phase space 
coordinates of all the particles that have been traced over. 
Similarly, also rj in Eq. ( 10 1 depends on variables associ- 



ated with particles that are being traced over. Therefore, 
also those contributions have to be recorded explicitly 
during the removal of particles from the system. 

In order to illustrate the present method we will apply 
it to the investigation of decoherence rates of superpo- 
sitions of harmonic oscillator coherent states \a) + \(3). 
The decoherence rate is predicted [TUJ [11] and experi- 
mentally verified |12l 113] to grow quadratically with in- 
creasing phase-space separation 5x between |a) and \0). 
However, this quadratic dependence holds only for small 
phase space separations, and once the separation is much 
larger then the range of the interaction potential a sat- 
uration is expected [HJ [TS]. We consider the full 3- 
dimcnsional problem with a thermal environment of mu- 
tually noniteracting particles that, however, interact with 
the oscillator via a short range Gaussian interaction. The 
environment is dilute so that only two-body interactions 
are taken into account. A thermal environment is real- 
ized via an average over 2000 realizations of an environ- 
ment consisting of 1500 particles. Due to the possibility 
of adding and removing particles this requires the inte- 
gration of 182 coupled differential equations, whereas the 
simulation of the entire system with all particles present 
during the entire integration would yield 8 x 10 7 differ- 
ential equations. 



FIG. 4: Decay constants for the coherence as in Fig. [3] as 
function of the initial relative phase space separation <5x be- 
tween | a) and \/3). For small 8x he decay constants grow 
quadratically with the phase space separation. Starting with 
8z exceeding the interaction range (25 oscillator length), a de- 
viation from the quadratic growth (dashed line) sets in, and 
the behavior is described better by a Gaussian (solid line). 



Fig. [3] shows the decay of the Hilbert Schmidt norm 
Q a p 5 the magnitude of coherence as function of 
time for different initial phase space separations 10, 20, 
30 and 40 oscillator lengths. The individual curves show 
an exponential decay exp(— jt) (with time given in mul- 
tiples of the oscillator period) with some remaining noise 
due to the average over a finite sample of environment re- 
alizations, and allow the extraction of a decay constant 7. 
Fig |4] displays this decay constant as function of the ini- 
tial phase space separation. Even for <5x vanishing, there 
is some decoherence with a decay constant of about 10 -5 
since a coherent state is not an eigenstate to the interac- 
tion potential, but the decay time for this case is signif- 
icantly longer than that for a coherent superposition of 
two coherent states. For small values of <5x, we recover 
the predicted quadratic increase of the decay constant, 
and once the phase space separation exceeds the range 
of the interaction potential (25 oscillator lengths in this 
case), the increase with <5x gets slower and, as expected, 
saturation sets in. 

The present techniques allow treatment of general situ- 
ations of quantum systems embedded in an environment. 
The range of applicability is basically limited only by the 
system-potentials, which should vary slowly over the typ- 
ical width of the wave packets that are used. Since wave 
packets spread spatially over time in many systems, there 
comes a time after which the present approximations be- 
come questionable. Nevertheless, since decoherence typ- 
ically takes place on time-scales that are significantly 
shorter than other damping phenomena, it is rather the 
short-term than long-term behavior that is of importance 
for open quantum systems. It is also important to note 
in this context that an environment interaction can re- 
sult in a narrowing of the system wave packet, so that 
with an environment there is an effect counteracting the 
spreading. 
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